Overview

Make three sets of plots for Around the island fDOM presentation

Indices: HIX, A, T, B

Details: Make colors pop (green to red) and drop outliers to get nice spread.

  1. Map

  2. Correlations to ‘distance to shore’

  3. Correlations to silicate, N+N, Phosphorous

library(ggplot2)
library(data.table)


theme_set(theme_minimal())

l_dat = fread("data/clean/ati_fdom_long.csv")

dat = fread("data/clean/ati_fdom_clean.csv")

meta = fread("data/sample_data/islandwide_sites_allMetadata.csv")


inds = c(
  'CobleA',
  'CobleB',
  'CobleT',
  'HIX'
)

# log transform all indices
l_dat$log_index_value = log(l_dat$index_value)

# set factor levels for Habitat

Summary of samples collected

site_meta = meta[!is.na(as.numeric(Site)) & May2021 == "Yes"]
## Warning in eval(.massagei(isub), x, ienv): NAs introduced by coercion
site_meta[!(Site %in% dat$Site) & May2021 == "Yes"]
##     Site Jan2016 May2016 July2016 July2019 Aug2020 May2021  Latitude Longitude
##  1:   30     Yes     Yes      Yes     <NA>     Yes     Yes -17.59838 -149.8108
##  2:   97     Yes     Yes      Yes     <NA>     Yes     Yes -17.54098 -149.9068
##  3:  122     Yes     Yes      Yes      Yes     Yes     Yes -17.48485 -149.8962
##  4:  125     Yes     Yes      Yes     <NA>     Yes     Yes -17.49020 -149.8827
##  5:  129     Yes     Yes      Yes     <NA>     Yes     Yes -17.48973 -149.8727
##  6:  141     Yes     Yes     <NA>     <NA>     Yes     Yes -17.49316 -149.8674
##  7:  156     Yes     Yes     <NA>     <NA>     Yes     Yes -17.48503 -149.8313
##  8:  207    <NA>    <NA>     <NA>     <NA>     Yes     Yes        NA        NA
##  9:  208    <NA>    <NA>     <NA>     <NA>     Yes     Yes        NA        NA
## 10:   15    <NA>    <NA>     <NA>     <NA>    <NA>     Yes        NA        NA
## 11:  210    <NA>    <NA>     <NA>     <NA>    <NA>     Yes        NA        NA
##     Distsance.to.crest Distance_to_shore Distance_to_pass
##  1:             716.16            718.52          3507.52
##  2:             194.05            563.28           160.22
##  3:             162.06            576.72            88.94
##  4:             789.48            137.46          1347.58
##  5:             618.38            232.22           930.83
##  6:             729.33            123.28           589.73
##  7:             843.02             44.33           787.09
##  8:                 NA                NA               NA
##  9:                 NA                NA               NA
## 10:                 NA                NA               NA
## 11:                 NA                NA               NA
##     Distance_to_deep_lagoon_water Distance_to_population_center
##  1:                        389.04                        718.52
##  2:                         99.38                       7334.53
##  3:                        274.81                       2147.15
##  4:                        371.65                        592.26
##  5:                         65.99                        232.22
##  6:                         31.14                        123.28
##  7:                         19.34                       2786.91
##  8:                            NA                            NA
##  9:                            NA                            NA
## 10:                            NA                            NA
## 11:                            NA                            NA
##     Population_center Island_shore       Habitat   Lagoon
##  1:           Ma'atea           SW    Mid lagoon   Maatea
##  2:    Vai'anae/Atiha           SW     Reef pass  Haapiti
##  3:         Papeto'ai            N     Reef pass    Hauru
##  4:         Papeto'ai            N Fringing reef Papetoai
##  5:         Papeto'ai            N Fringing reef Papetoai
##  6:         Papeto'ai            N Fringing reef Papetoai
##  7:            Paopao            N Fringing reef  Pihaena
##  8:              <NA>         <NA>          <NA>     <NA>
##  9:              <NA>         <NA>          <NA>     <NA>
## 10:              <NA>         <NA>          <NA>     <NA>
## 11:              <NA>         <NA>          <NA>     <NA>

plots

Map

library(sf)
## Linking to GEOS 3.10.2, GDAL 3.4.3, PROJ 8.2.1; sf_use_s2() is TRUE
library(sp)
library(OpenStreetMap)
library(viridisLite)

map_long_dat = l_dat[!is.na(Latitude)]

# make spacial points data frame
map_sp = SpatialPointsDataFrame(data = map_long_dat,
                       coords = list(map_long_dat$Latitude,
                                  map_long_dat$Longitude))


# Check geographic range of sampling points
limits = c(
  min(map_long_dat$Longitude),
  min(map_long_dat$Latitude), 
  max(map_long_dat$Longitude),
  max(map_long_dat$Latitude) 
)

# define a bounding box with a small cushion around the minimum and maximum
bbox = list(
  xmin = limits[1] - 0.03,
  ymin = limits[2] - 0.04,
  xmax = limits[3] + 0.03,
  ymax = limits[4] + 0.04
)

# get basemap
sa_map <- openmap(c(bbox$ymax, bbox$xmin),
                  c(bbox$ymin, bbox$xmax),
                  type = "stamen-terrain",
                  mergeTiles = TRUE)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition
sa_map2 <- openproj(sa_map)


mo_map = function(an_index, outlier_n){
  
  # remove highest and lowest n samples
  ind_map_dat = map_long_dat[index_name == an_index &
                               !(is.infinite(log_index_value))][order(log_index_value)]
  ind_map_dat = head(ind_map_dat, n = -(outlier_n))
  ind_map_dat = tail(ind_map_dat, n = -(outlier_n))

  sa_map2_plt = OpenStreetMap::autoplot.OpenStreetMap(sa_map2)+
    geom_point(data = ind_map_dat,
               aes(x = Longitude,
                   y = Latitude,
                   color = log_index_value))+
    labs(title = paste("Log",an_index),
         x = "Lon",
         y = "Lat")+
    scale_color_gradient(low = "green", high = "red")+
    guides(fill=guide_legend(title=NULL))

print(sa_map2_plt)

ggsave(filename = paste0("plots/exploration/ATI/map_",an_index,".png"), plot = sa_map2_plt)
  
}
 

lapply(inds, mo_map, outlier_n = 10)
## Saving 10 x 5 in image

## Saving 10 x 5 in image

## Saving 10 x 5 in image

## Saving 10 x 5 in image

## [[1]]
## [1] "plots/exploration/ATI/map_CobleA.png"
## 
## [[2]]
## [1] "plots/exploration/ATI/map_CobleB.png"
## 
## [[3]]
## [1] "plots/exploration/ATI/map_CobleT.png"
## 
## [[4]]
## [1] "plots/exploration/ATI/map_HIX.png"

Distance to Shore

# function for making correlation plots
cor_plot = function(ind,
                    cor_var,
                    outlier_n,
                    cor_var_outlier_n){
  
  #subset and drop outliers
  ind_dat = l_dat[index_name == ind & !is.infinite(log_index_value)][order(index_value)]
  ind_dat = head(ind_dat, n = -(outlier_n))
  
  if(!is.null(cor_var_outlier_n)){
  ind_dat = head(ind_dat[order(ind_dat[[cor_var]])], n = -(cor_var_outlier_n))
  
  }
  
  # model
  ind_lm = lm(as.formula(paste("log_index_value ~", cor_var)),
              data = ind_dat)
  ind_sum = summary(ind_lm)
  
  p = ggplot(data = ind_dat, aes_string(x = cor_var,
                                y = "log_index_value"))+
    geom_point(aes(color = Habitat))+
    
    geom_abline(slope = ind_lm$coefficients[2],
                intercept = ind_lm$coefficients[1])+
    
    labs(y = paste( "Log",ind),
         subtitle = paste("R sq. =",
                          round(ind_sum$r.squared, 3)))

  ggsave(paste0("plots/exploration/ATI/",cor_var,"_",ind,".png"))
  
  return(p)
  
}
# apply correlation plots for distance to shore
d_shore_ps = lapply(inds,
                    cor_plot,
                    cor_var = "Distance_to_shore",
                    outlier_n = 10,
                    cor_var_outlier_n = NULL)
## Saving 10 x 5 in image
## Warning: Removed 7 rows containing missing values (geom_point).
## Saving 10 x 5 in image
## Warning: Removed 7 rows containing missing values (geom_point).
## Saving 10 x 5 in image
## Warning: Removed 8 rows containing missing values (geom_point).
## Saving 10 x 5 in image
## Warning: Removed 8 rows containing missing values (geom_point).
d_shore_ps
## [[1]]
## Warning: Removed 7 rows containing missing values (geom_point).

## 
## [[2]]
## Warning: Removed 7 rows containing missing values (geom_point).

## 
## [[3]]
## Warning: Removed 8 rows containing missing values (geom_point).

## 
## [[4]]
## Warning: Removed 8 rows containing missing values (geom_point).

Nutrients

# apply correlation to nutrients

# N+N
NplusN_ps = lapply(inds,
                    cor_plot,
                    cor_var = "Nitrite_plus_Nitrate",
                    outlier_n = 10,
                    cor_var_outlier_n = 2)
## Saving 10 x 5 in image
## Saving 10 x 5 in image
## Saving 10 x 5 in image
## Saving 10 x 5 in image
NplusN_ps
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

# Phosphate
Phos_ps = lapply(inds,
                    cor_plot,
                    cor_var = "Phosphate",
                    outlier_n = 10,
                    cor_var_outlier_n = 2)
## Saving 10 x 5 in image
## Saving 10 x 5 in image
## Saving 10 x 5 in image
## Saving 10 x 5 in image
Phos_ps
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

# Silicate
Sil_ps = lapply(inds,
                    cor_plot,
                    cor_var = "Silicate",
                    outlier_n = 10,
                    cor_var_outlier_n = 4)
## Saving 10 x 5 in image
## Saving 10 x 5 in image
## Saving 10 x 5 in image
## Saving 10 x 5 in image
Sil_ps
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

Habitat Boxplots

hab_box = function(ind){
  
p = ggplot(l_dat[index_name == ind], aes(x = index_name,
                                         y = log_index_value,
                                         color = factor(Habitat, levels = c("Reef pass",
                                                                            "Reef crest",
                                                                            "Mid lagoon",
                                                                            "Bay",
                                                                            "Fringing reef",
                                                                            ""))))+
      geom_boxplot()+
      labs(title = paste0("log ",ind))+
      facet_wrap(vars(Island_shore), nrow = 2)+
      theme(legend.title=element_blank())


print(p)
  
ggsave(paste0("plots/exploration/ATI/Habitat_boxplot_",ind,".png"))
  
}


lapply(inds, hab_box)
## Saving 10 x 5 in image

## Saving 10 x 5 in image

## Warning: Removed 16 rows containing non-finite values (stat_boxplot).
## Saving 10 x 5 in image
## Warning: Removed 16 rows containing non-finite values (stat_boxplot).

## Warning: Removed 1 rows containing non-finite values (stat_boxplot).
## Saving 10 x 5 in image
## Warning: Removed 1 rows containing non-finite values (stat_boxplot).

## [[1]]
## [1] "plots/exploration/ATI/Habitat_boxplot_CobleA.png"
## 
## [[2]]
## [1] "plots/exploration/ATI/Habitat_boxplot_CobleB.png"
## 
## [[3]]
## [1] "plots/exploration/ATI/Habitat_boxplot_CobleT.png"
## 
## [[4]]
## [1] "plots/exploration/ATI/Habitat_boxplot_HIX.png"